Hyperbolic chaos at blinking coupling of noisy oscillators 
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We study an ensemble of identical noisy phase oscillators with a blinking mean-field coupling, 
where one-cluster and two-cluster synchronous states alternate. In the thermodynamic limit the 
population is described by a nonhnear Fokker-Planck equation. We show that the dynamics of 
the order parameters demonstrates hyperbohc chaos. The chaoticity manifests itself in phases of 
the complex mean field, which obey a strongly chaotic Bernoulli map. Hyperbolicity is confirmed 
by numerical tests based on the calculations of relevant invariant Lyapunov vectors and Lyapunov 
exponents. We show how the chaotic dynamics of the phases is slightly smeared by finite-size 
fiuctuations. 
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I. INTRODUCTION 

Dynamics of large populations of coupled oscillators 
attracted large interest recently. Problems of this type 
appear in studies of Josephson junctions, lasers, ensem- 
bles of neurons, cell populations, and many other fields. 
From a more general perspective, studies of such a sys- 
tem allows one to shed light on a long-standing problem 
of an interrelation between microscopic and macroscopic 
dynamics. Indeed, the mostly studied nontrivial effect 
in the ensembles of globally coupled oscillators is their 
synchronization, which can be considered as a nonequi- 
librium phase transition [1. 2]. Remarkably, in some sit- 
uations one can explicitly derive the dynamics of global 
variables (order parameters), in terms of which the syn- 
chronization transition is a bifurcation from zero equilib- 
rium to a non-trivial state [3]. 

Different interrelations between regularity properties 
on micro- and macro-levels (i.e. for individual oscillators 
and for mean fields) have been reported in the literature. 
For example, chaotic micro-oscillators being coupled may 
lead to periodic mean fields [4-8]. On the other hand, 
coupled periodic oscillators may produce chaotic mean 
fields [9-11]. While description of populations of deter- 
ministic oscillators remains a challenging task, there ex- 
ists a nice framework for ensembles of noise-driven oscilla- 
tors. The behavior in the thermodynamic limit can be de- 
scribed self-consistently with a nonlinear Fokker-Planck 
equation, first suggested by Desai and Zwanzig ]12] in 
the context of globally coupled noisy bistable oscillators 
(see also [13[). This approach has been then applied to 
noisy periodic oscillators in Refs. [14, 15]. 

In this paper we consider a population of noisy oscilla- 
tors subject to blinking, time-periodic coupling ]16-18]. 
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This is a minimal generalization of the simplest model 
with constant coupling, which demonstrates simple syn- 
chronization patterns only. We will show that with a 
blinking coupling, where on different periods of the to- 
tal cycle different synchronous modes emerge, the total 
dynamics demonstrates highly nontrivial regime of phase 
hyperbolic chaos. In this regime the phase of the complex 
order parameter obeys a doubling Bernoulli map, which 
has strong chaotic properties, and is, contrary to many 
other models of chaos, structurally stable in respect to 
perturbations. 

Our treatment of the dynamics of the complex mean 
fields follows recent studies [19, 20], where a general 
framework for mechanisms of hyperbolic chaos in coupled 
oscillators has been developed. In the simplest setup, 
two oscillators that are excited alternately can interact 
in a way to influence the phases of each others at the 
stages where these oscillators pass through a Hopf bifur- 
cation. More precisely, one needs that the appearing field 
is forced by a higher harmonic of another oscillator, then 
a transformation of the phase — > ncf) with |n| > 2 can 
occur. As the result the phase obeys a Bernoulli-type, 
uniformly expanding map, and the whole strange attrac- 
tor is of famous Smale- Williams type [21]. The described 
mechanism of phase multiplication generating hyperbolic 
chaos is quite generic. The straightforward way of its 
implementation is to consider two non-autonomous self 
oscillators with different natural frequencies and appro- 
priate coupling [22]. Further analysis revealed a possi- 
bility to identify this mechanism in a system of two or 
three coupled autonomous oscillators [23]. Recently it 
was described, how to observe the hyperbolic chaos in a 
spatially extended system, as a result of interaction of 
the Turing modes [24]. In paper ]11] an alternation be- 
tween synchronized and desynchronized regimes in two 
ensembles of non-identical oscillators subjected to the 
Kuramoto transition was shown to possess a collective 
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hyperbolic chaos of complex order parameters. 

In this paper we study an ensemble of identical stochas- 
tic phase oscillators coupled via the mean fields. We con- 
sider a situation of blinking coupling, where different syn- 
chronization patterns, one with one cluster and another 
with two clusters, alternate. We demonstrate, both bas- 
ing on the nonlinear Fokker-Planck equation, and on the 
direct simulation of a large population that the phases 
of these patterns demonstrate hyperbolic chaos. More- 
over, we study the Lyapunov exponents of the system 
and apply numerical criteria based on Lyapunov vectors 
to verify hyperbolicity. 



II. ENSEMBLE OF COUPLED NOISY 
OSCILLATORS 

A. Microscopic equations 

Our basic model is an ensemble of K globally coupled 
identical limit-cycle oscillators with additive noise. In 
the phase description the ensemble is governed by the 
following set of equations: 

1 ^ 

0fe = fiM + ]f E - M + 'rUt), (1) 

where k = 1, . . . , K, (pk is the phase of a fc-th oscillator 
(taken in the reference frame rotating with the basic fre- 
quency), and the independent random variables ^fc(t) de- 
scribe a white Gaussian noise: {(,k{t)^j {t + r)} = (5jfc(5(r). 
Functions and /(0) (both are 27r-periodic) describe, 
respectively, effects of global coupling and of the exter- 
nal periodic forcing at the oscillator frequency and/or its 
harmonics. It is convenient to define these functions by 
means of the Fourier decompositions 

oo oo 
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Since the functions / and q are real, i^_„ = F*, and 
Q~n = Qn^ where asterisk denotes complex conjugation. 

Given a state of the ensemble, one can determine the 
following mean fields playing roles of order parameters: 

A,(^) = (e'^*) = -if;e'^^^ (3) 

k=l 

where £ = 1,2,..., and = 1. Notice that A-i = A^. 
In what follows we shall denote a,rg Ai = 



B. Nonlinear Fokker-Plank equation 

In the thermodynamic limit K ^ oo the ensemble can 
be described by the density of the probability distribution 



v{(f), t). The mean fields in this case are expressed as 

Mt)^ r ^{<l>,t)e''Ucl). (4) 
Jo 

The probability density can be decomposed as Fourier 
series: 

oo 

«(<^'^) = ^ E Mt)e-''t (5) 

e=-oo 

i.e.. the order parameters A^(t) are just coefficients of the 
Fourier modes for the probability density. 

The dynamics of the probability density can be de- 
scribed by the nonlinear Fokker-Planck equation [12, 13]. 
For the system (1), accounting the expressions (2), we can 
write out this equation as 

n— — oo 

Now we substitute the Fourier decomposition (5) for 
v{(j), t), then multiply the equation by e"'^"^ and integrate 
it over 2tt. The resulting equations for Ag read: 

Ai = ie 22 {F-n+QnAn)Ae-n- ^Ag. (7) 

n— — oo 

C. Elementary synchronization dynamics 

The set of equations (7) allows one a simple qualitative 
description of basic synchronization phenomena in the 
ensemble. If only the global coupling is present (F„ = 0), 
and q{(j)) has one Fourier component Qn — — ^k„, then 
the non-synchronized state Aj = 6j_o remains stable until 
the coupling passes the synchronization threshold, 

Kn > na^. (8) 

Above this threshold, the mode A„ becomes unstable, 
and a stationary solution of (7) containing modes Ain 
with £ = 1,2,3,... establishes. This corresponds to the 
appearance of an n-clustered state of the ensemble, where 
the probability density has maxima at cj)^^"^ + fc27r/n, 
and i/)^*'-' is an arbitrary phase shift determined by ini- 
tial conditions. In contrast, without the global coupling 
(Qn ==0), if only the external forcing is present contain- 
ing one Fourier component {Fn = — ^e), then for small 
e an n-cluster state with A^ ~ — en~^o'~'^ appears. The 
phase of this cluster is fixed, been determined by the ap- 
plied forcing. 

III. CLUSTER EXCHANGE DYNAMICS AT 
BLINKING COUPLING 

The main goal of this paper is to describe macroscopic 
chaos appearing when the coupling in the ensemble is 
blinking (cf. [16]). 
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In this section we consider the effect of bhnking cou- 
pling qualitatively, basing on the elementary synchro- 
nization dynamics described in the previous section. 
First, we make a particular choice for the global coupling 
terms, assuming alternating couplings Qi and Q2'. 



stage 2 to stage 1, we obtain the Bernoulli- type map 



l'«l,2 



n:2 



(9) 



where 



(ki,K2) 



(ki,0), \imT <t <mT + T/2, 
(0,K2), ilmT + T/2 <t <mT + T. 



(10) 

Here T is the total period of blinking; it is assumed to 
be large enough compared to characteristic time of clus- 
ter formation or decay. In the first half of this period 
(we call it stage 1) only Qi is present, correspondingly a 
one-cluster state develops; during the second half of the 
period (stage 2) only Q2 is present, so in this stage the 
two-cluster state develops (see Fig. 1 below). 

For the following consideration, the phases of the clus- 
ters are important. As the one-cluster state contains 
all modes A„, n = 1,2,3,..., at the beginning of the 
stage 2 the amplitude A2 ~ A\ is finite. On the stage 2, 
due to the appropriate coupling the instability appears 
associated with formation of the two-cluster state, and 
the initial phase $2 = 2<I>i occurs to be definite, deter- 
mined by the preceding evolution. Thus, the two-cluster 
state accepts this phase from the former one-cluster state. 
At the end of stage 2, only modes with even index A„, 
n = 2, 4, 6, . . . are present (as T is large enough, the odd 
modes decay effectively during the stage 2). Thus, at the 
beginning of the next stage 1 there is no initial amplitude 
Ai, and its growth due to the coupling Qi would start 
only from random initial fluctuations, and the appearing 
one-cluster state would have a random phase. 

The situation changes if a small regular external forc- 
ing is present. As one can see from (7), the forcing Fn 
provides appearance of terms iF^nA^-n, thus produc- 
ing combinational modes. 

Let us consider the effect of the mode F^ — At 
the stage 1 it produces combinational modes with all in- 
dices, but because these modes are already present in 
the distribution, the effect of small e is negligible. At 
the stage 2, however, the mode Ai is produced by terms 
^3*^2 +-^3^4 in (7). Because A2 is typically much larger 
than An, the result of this interaction is the appearance 
of a small but notable amount of mode 1 with Ai ^ tA^ 
at the end of the stage 2. 

This circumstance changes essentially the starting con- 
ditions for the evolution on the stage 1 compared to the 
case e = 0: now the growing mode Ai develops from the 
seed ~ 6^2, and so its phase will be $1 = — <I>2- Now 
the phases of cluster states arising at all stages of evolu- 
tion are non-random, but depend on the previous phases 
in a deterministic way. Combining the transformations 
<f>2 = 2(f>i from stage 1 to stage 2, and $1 = — $2 from 



<i>i(TO -I- 1) 



-2<I>i(m) 



(11) 



at successive periods of the blinking coupling. 

The map (11) for the phase transformation is uniformly 
expanding and hyperbolic, and we conclude from these 
qualitative arguments that the cluster patterns for blink- 
ing coupling of noisy oscillators in the situation we have 
considered will demonstrate hyperbolic chaos in the dy- 
namics of the order parameters (the modes of the distri- 
bution density). 



IV. ILLUSTRATIONS AND 
CHARACTERIZATION OF THE DYNAMICS 

In this section we confirm numerically the existence of 
the hyperbolic chaos in the setup described above. The 
mode equations take the form: 



A, 



^ ^ K„ [AfiAi-n A—nAi-f-n) 



i+3 



(12) 



The corresponding Fokker-Plank equation reads 
dv d 



dt 




(13) 



where (sin[n(?/; — 0)])^ = J^^ ^("0, t) sm[n('ip—(f))]d'ip. Ac- 
counting this form, we can reconstruct the respective mi- 
croscopic equation (1) for the ensemble of phase oscilla- 
tors as 

2 , K 

= esin(30fc) + ^ ^ sin[n((/)j - (t>k)] + (Jik{t). 

n—1 j—^ 

(14) 

In all these equations the blinking coupling K„(t) is given 
by expression (10). 

The essential control parameters both for the ensem- 
ble (14). and for the Fokker-Plank system (12) are ki and 
K2, which have to satisfy the condition (8) for clusters to 
grow at the corresponding stages. 



A. Nonlinear Fokker-Plank equation 

First, we present numerics for the Fokker-Plank sys- 
tem (12). Since this system appears as a spectral de- 
composition of the nonlinear integro-differential equa- 
tion (13), it is stiff, and we integrate it numerically using 
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the method of exponential time differencing, see Ref. [25]. 
As fohows from Eq. (3), \Ai\ < 1. Due to the terms pro- 
portional to £~'^, the amplitudes of Ae decay fast with 
£. Direct verification shows that already for i = 5 the 
magnitude is small enough, jAsj ^ 1. Thus we cut the 
infinite set of Eqs. (12) at ^cut = 10. 

Figure 1(a) illustrates the temporal behavior of the 
probability density «(</), t) computed according to Eq. (5). 
At the stage 1 (0 < t < 75) the 1-cluster develops. Then, 
after the switch to the stage 2 at i = T/2 = 75, the 
density acquires two humps. One part of the 2-cluster is 
located just at the site of the former 1-cluster, and the 
second one appears shifted by tt. As we already discussed, 
this situation corresponds to the quadratic relation for 
the complex amplitudes A2 ~ A^ and to the doubling 
of the argument. At the initial epoch of the next stage 

1 at t > 150; the 2-cluster disappears, and the 1-cluster 
reappears. Notice that the hump emerges now at a site 
which is different both from that of the previous 2-cluster, 
and that of the former 1-cluster. It is so because of the 
presence of a small force f{4>) corresponding to the term 
proportional to e in Eqs. (12). 

Figures 1(b) and (c) show magnitudes and arguments 
of An- At stage 1 all the magnitudes \An\ grow and 
saturate at some level; in this state the main mode Ai 
dominates. At stage 2 the mode A2 grows, together with 
its harmonics, while all odd modes decay. 

Let us now focus on the phase dynamics. The phase 
<i>2 does not change at the transition from stage 1 to stage 

2 at t = 75. At the beginning of stage 2 the phase of the 
first mode $1 is nearly constant while the amplitude of 
this mode remains relatively large; however at i « 100 the 
amplitude drops to level ^ e, and now this mode becomes 
to be driven by eAj. Here, the mode Ai accepts the phase 
—$2 that can be seen in Fig. 1(c) as a sharp transition 
around t fa 110. After this event, the phase of the mode 
Ai does not vary significantly. At t = T = 150 the new 
stage 1 starts, here the amplitude of the previously active 
mode A2 decays rapidly, but it reappears soon as the 
second harmonics of the mode 1, accepting its doubled 
phase 2$i (see the jump of $2 at t « 155 followed by the 
growth of 1^2!). 

Although the dynamics of the phases and the ampli- 
tudes observed in Fig. l(a,b,c) nicely corresponds to the 
qualitative picture of section II C, it is important to ver- 
ify that the phase evolution really follows the Bernoulli- 
like map (11). In Fig. 2(a) a diagram is shown for the 
values of argument of the complex variable Ai recorded 
stroboscopically at t = t„ = (n + ■^)T. Indeed, the 
plot demonstrates the doubling of the phase, as expected. 
Thus, actually, the phases of the alternately arising clus- 
ters evolve according to the expanding circle map that 
implies the hyperbolic chaos. 



B. Ensemble of noisy oscillators 

Here we report on the results of a straightforward sim- 
ulation of the ensemble dynamics according to Eqs. (14). 
We perform numerical solutions for these stochastic dif- 
ferential equations employing a special version of the 
Runge-Kutta-type algorithm suggested in [26]. Fig- 
ure 1(d) illustrates the temporal behavior of the prob- 
ability density v{(j),t) computed from the instantaneous 
distributions of the phases 0^ for the elements of the 
ensemble, and the panels (e) and (f) demonstrate the 
dynamics of the mean-field variables A„(<). The param- 
eters are the same as in the previous simulations of the 
Fokker-Plank equation. Except for the obvious difference 
in the initial conditions, the distributions v{(j),t) and the 
magnitudes of A„(t) at their high levels behave visually 
similar. The dynamics of the phases is definitely more 
noisy than that in Fig. 1(c), especially in the regions 
where the magnitudes of the corresponding modes are 
small: in these regions even small finite-size fluctuations 
produce a large effect on the phases. When the magni- 
tude of ^i(t) is high (see for example the time interval 
30 < t < 75); its argument does not respond to the 
noise. But when Ai gets small, the argument fluctuates 
(see the time interval 110 < i < 150). The same is true 
for A2. The most harmful are the fluctuations of $1 be- 
cause this variable inherits the argument value from A2 
when I All is small. If the fluctuations are strong enough, 
they wash out the correct value of the argument, and the 
phase multiplication mechanism ceases to operate prop- 
erly. The fluctuations of A2 are not so essential since 
they occur only in epochs, when the component A2 does 
not play an important role. One can see from the fig- 
ure that the phase transfer happens when A2 has a large 
magnitude and thus is not affected by the noise (see the 
time interval 180 < t < 225). 

Figure 2(b) shows the diagram for the stroboscopically 
recorded values of $1 corresponding to the dynamical be- 
havior illustrated in Fig. l(d,e,f). One can observe that 
the phases demonstrate the expected Bernoulli-like dou- 
bling mapping in average, though the fluctuations pro- 
duce a noticeable statistical widening (which decreases 
as the size of the ensemble K grows). 

C. Lyapunov exponent and hyperbolicity test 

Now we return to the Fokker-Plank system (12) and 
discuss how the Lyapunov exponents (LEs) depend on ki 
and K2- Also we test hyperbolicity of the chaotic dynam- 
ics of the system (12). 

We calculate the Lyapunov exponents in the determin- 
istic system (12) with the standard Benettin algorithm. 
Because the system is periodically driven, we report as 
LEs the values A = AT, i.e. the "Lyapunov multipli- 
ers" over the period. In all cases we have found at most 
one positive LE, the others were negative. Figure 3(a) 
shows the largest LE versus parameter ki; in a large 
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Figure 1. (color online) a) Probability density v{(fi,t) determined by (5), where An{t) are computed according to Eq. (12); d) 
distribution for the ensemble of c/tk, see Eq. (14), where the total number of the oscillators is K — 10^. Darker colors correspond 
to higher values of v[<j),t). Five periods of switching of Ki,2(t) are shown, b) and c) Absolute values and arguments of An 
computed according to Eq. (12); e) and f) the same computed via Eq. (3) for the ensemble. \An\ are plotted in logarithmic 
scale. Parameters: ki = 1.3, k.2 = 2.3, T = 150, e = 0.01, cr = 1, Ic^it = 10. 




Figure 2. a) The plot tpn+i vs. ipn where ip„ = axg A\{nT + 
^T), and Ai is a solution of Eqs. (12). b) The same diagram 
for the ensemble (14), where A\ is computed according to 
Eq. (3). Parameters correspond to Fig. 1. 



range of this parameter it is positive and close to the 
expected value ln2 0.69314. The other exponents 
are all negative, and their magnitudes are much higher. 
For example, at ki = 1.3, and K2 = 2.3 we have = 
{0.688, -17.355, -197.507, -245.781, -705.748}, and the 
corresponding Kaplan- Yorke dimension is 1.04. 

The positive Lyapunov exponent varies slowly depend- 
ing on the parameter, but no tips or dips, like in many 
non-hyperbolic systems, are observed. This may be re- 
garded as a manifestation of the structural stability in- 
trinsic to the hyperbolic dynamics. Figure 4(a) repre- 
sents dependence of A on K2, while ki is constant. Again, 



the same features are observed: the largest Lyapunov ex- 
ponent is close to the value In 2, and it does not demon- 
strate any notable variation within a wide range of the 
parameter. 

To test hyperbolicity, we employed the method which 
has been specially developed for high-dimensional sys- 
tems, see Ref. [27] for details and Ref. [28] for the math- 
ematical background. In brief, the method includes the 
computation of the first k orthogonalized Lyapunov vec- 
tors moving forward in time (so called Gram-Schmidt 
vectors or backward Lyapunov vectors [28]), and also the 
first k orthogonalized Lyapunov vectors obtained moving 
backward in time from the conjugate variational equa- 
tions (forward Lyapunov vectors) ; here k is the dimension 
of the unstable manifold. Given these vectors, a. k x k 
matrix of their scalar products is built. Its smallest sin- 
gular value dfc is cosine of the angle between the expand- 
ing tangent subspace and the orthogonal conrplemeirt to 
the contracting tangent subspace, and thus dk can be 
an indicator of hyperbolicity. For a non-hyperbolic case, 
the dk values vanish somewhere along the trajectories 
(as the tangencies between the associated expanding and 
coirtracting subspaces occur), but for situations of the hy- 
perbolic dynamics the distribution of dk is well separated 
from zero. 

In Fig. 3(b) we illustrate the results of this hyperbol- 
icity test. In our case the expanding subspace is one- 
dimensional. Thus, for each ki we compute 10^ points on 
the attractor of the stroboscopic map, and evaluate the 
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D. Lyapunov exponents for the micro-level 
dynamics 

It is instructive to compare the macroscopic LEs for the 
dynamics of the order parameters with microscopic LEs 
describing stabihty of the dynamics of individual oscilla- 
tors. A direct approach to computation of the full spec- 
trum of Lyapunov exponents for the ensemble of phase 
oscillators (14) can easily exhaust a computer memory 
and take an extremely long time. However, there is a 
simple way to estimate microscopic LEs assuming their 
decoupling from the macroscopic dynamics. Neglecting 
small e in Eq. (3), we can rewrite this equation as 



Figure 3. (a) The first Lyapunov exponent, and (b) the 
lower and upper boundaries of the distribution of di vs. m. 
The data relate to the stroboscopic map corresponding to the 
Fokker-Plank system (12) at tn = {n + ■^)T . Distributions 
of di for each ki were computed for 10* attractor points. 
K2 = 2.3, other parameters correspond to Fig. 1. 



I ^ Ai„(A„e-™'^'= - A:e'"^-) + aCkit) • (15) 



Let us now assume that on stages 1 and 2 the order pa- 
rameters Ai and A2 are roughly constant. Then the vari- 
ations of the phases satisfy 
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Figure 4. The same quantities as in Fig. 3 against K2, while 
Kl = 1.3. 



k dt 



a: 



(16) 



Now the microscopic LE is evaluated as an average of 
the r.h.s. of Eq. (16) multiplied by T/2, which yields 

Amic ~ ^1" '^'^nl^nP- Substituting here the ap- 

proximate expressions for stationary magnitudes of the 
order parameters, valid for small subcritical couplings 
Ki > cr^, K2 > 2(T^, we obtain the following estimate: 



A, 



Kl 



K2 - 2a' 

[— 

K2 



(17) 



For the parameters used above, a ~ 1,ki ~ 1.3, K2 = 
2.3, T = 150, this yields A,nic « —113 in a reasonable 
agreement with the numerics. 

Thus, contrary to many situations where chaos oc- 
curs in the micro-level description, while the macro- 
parameters manifest rather regular behavior, here we ob- 
serve an impressively opposite case: the micro-dynamics 
is stable (because Amic < 0), but at the macro-level the 
hyperbolic chaos takes place. 



indicator di at each of them. Then, the smallest and the 
largest values of di are plotted on the diagram, marked 
with crosses and pluses, respectively. The first observa- 
tion is that the smallest c?i is very far from zero; this 
confirms the hyperbolicity. Another noteworthy obser- 
vation is that the maximal and the minimal values of di 
are very close to each other. This means that the mutual 
orientations of the expanding and contracting directions 
remain practically unaltered on the entire attractor. This 
can be treated as another manifestation of the structural 
stability. Figure 4(b) shows the hyperbolicity indicator 
di in dependence on K2, it is again well separated from 
zero. 



V. CONCLUSION 

In this paper we have demonstrated that blinking cou- 
pling in populations of noisy oscillators may result in 
hyperbolic chaos of macroscopic mean fields. The blink- 
ing is between two simplest modes of coupling: one leads 
to a one-cluster synchronized state, another one leads to 
appearance of two clusters, shifted by tt. Correspond- 
ingly, at these couplings two different order parameters 
dominate: the usual Kuramoto order parameter for the 
one-cluster state, and the second-harmonics Daido order 
parameter for the two-cluster state. While the magni- 
tudes of these order parameters show no significant vari- 
ations from cycle to cycle, the phases obey a strongly 
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chaotic Bernoulli transformation, and thus demonstrate 
a hyperbolic chaos. Noteworthy, there is a realistic situ- 
ation where the two types of global coupling (one-cluster 
and two-cluster) can be observed: this is an ensemble of 
pendula hanging on a common beam. In such a configu- 
ration with two pendulum clocks Ch. Hyugens observed 
synchronization more than 300 years ago (see a transla- 
tion of his notes in [2]). Nowadays one reproduces these 
experiments with metronomes [29]. The horizontal mo- 
tions of the beam result in the one-cluster coupling, while 
the vertical motions lead to the two-cluster coupling [30] . 

Our analysis is mostly based on the consideration of 
the thermodynamic limit, where the population can be 
described by the nonlinear Fokker-Planck equation. We 
verify hyperbolic chaos in this integral-differential equa- 
tion, by direct simulations of equations for the modes and 
showing that their phases obey a Bernoulli map, and by 
checking that the stable and the unstable directions in 
the tangent space are never tangent to each other. For 
the original population a direct simulation gives a pic- 
ture very similar to that for the nonlinear Fokker-Planck 



equation, but the one-dimensional transformation of the 
phases looks like a noisy Bernoulli map, due to finite-size 
effects. 

In a more general perspective, the nonlinear Fokker- 
Planck equation, as an integro-differential equation with 
partial derivatives, is a representative of a class of deter- 
ministic distributed systems demonstrating pattern for- 
mation. Different cluster synchronization states corre- 
spond to different patterns; phases of clusters correspond 
to the spatial positions of the patterns. With this in- 
terpretation, our results demonstrate that an interaction 
between blinking patterns results in a chaotic relocation 
of their positions, moreover, this chaos is hyperbolic. A 
kind of such behavior was reported recently in a model 
of interacting Turing patterns with different wave num- 
bers [24]. The presented model based on the nonlinear 
Fokker-Plank equation provides another indication that 
the above mechanism of the hyperbolic chaos is realizable 
in quite generic circumstances. 
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